Regresar a la página principal
El análisis multivariante es conjunto de métodos estadísticos cuya finalidad es analizar simultáneamente conjuntos de datos multivariantes: hay varias variables medidas para cada caso.
Permiten un mejor entendimiento del fenómeno objeto de estudio, obteniendo información que los métodos univariantes y bivariantes son incapaces de conseguir.
Clasificación de las técnicas multivariantes
Las técnicas multivariantes se pueden clasificar según dos posibles criterios:
Métodos Dependientes. Se está interesado en la asociación entre las distintas variables, es decir, en las relaciones entre las mismas, donde parte de estas variables dependen o se miden en función de las otras. Subyace en ellos siempre un interés predictivo. Por ejemplo: regresión lineal multiple, análisis discriminante, modelos log-lineales, correlaciones canónicas, análisis multivariante de la varianza.
Métodos Independientes. Se está interesado en investigar las asociaciones que se presentan entre variables sin distinción de tipos entre ellas. Tienen un interés descriptivo. Por ejemplo: análisis factorial, escalado multidimensional, análisis de correspondencias, análisis cluster.
Cuando se recoge la información de una muestra de datos, lo más frecuente es tomar el mayor número posible de variables. Sin embargo, si tomamos demasiadas variables sobre un conjunto de objetos, por ejemplo 20 variables, tendremos que considerar \(\binom{20}{2}=180\) posibles coeficientes de correlación, si son 40 variables el número aumenta a 780. Por tanto, en este caso es difícil visualizar relaciones entre las variables.
Otro problema que se presenta es la fuerte correlación que muchas veces se presenta entre las variables: si tomamos demasiadas variables (cosa que en general sucede cuando no se sabe demasiado sobre los datos o sólo se tiene ánimo exploratorio), lo normal es que estén relacionadas o que midan lo mismo bajo distintos puntos de vista. Por ejemplo, en estudios médicos, la presión sanguínea a la salida del corazón y a la salida de los pulmones están fuertemente relacionadas.
Se hace necesario, pues, reducir el número de variables. Es importante resaltar el hecho de que el concepto de mayor información se relaciona con el de mayor variabilidad o varianza. Cuanto mayor sea la variabilidad de los datos (varianza) se considera que existe mayor información, lo cual está relacionado con el concepto de entropía.
Las nuevas variables son combinaciones lineales de las anteriores y se van construyendo según el orden de importancia en cuanto a la variabilidad total que recogen de la muestra. De modo ideal, se buscan \(p<m\) variables que sean combinaciones lineales de las \(m\) originales y que estén incorreladas, recogiendo la mayor parte de la información o variabilidad de los datos.
Matemáticamente
Supongamos que existe una muestra con \(n\) individuos para cada uno de los cuales se han medido \(m\) variables (aleatorias) \(F_j\). El PCA permite encontrar un número de factores subyacentes \(p < m\) que explican aproximadamente el valor de las \(m\) variables para cada individuo. El hecho de que existan estos \(p\) factores subyacentes puede interpretarse como una reducción de la dimensionalidad de los datos: donde antes necesitabamos \(m\) valores para caracterizar a cada individuo ahora nos bastan \(p\) valores. Cada uno de los \(p\) encontrados se llama componente principal.
Existen dos formas básicas de aplicar el ACP:
Método basado en la matriz de correlación, cuando los datos no son dimensionalmente homogéneos o el orden de magnitud de las variables aleatorias medidas no es el mismo.
Método basado en la matriz de covarianzas, que se usa cuando los datos son dimensionalmente homogéneos y presentan valores medios similares.
El primer paso para definir los componentes principales de p variables originales es encontrar una función lineal \(a_1^\prime y\), donde \(a_1\) es un vector de \(p\) constantes, para los vectores con varianza máxima.
\[ a_1^\prime y = a_{11}x_1 + a_{12}x_2 + \cdots + a_{1p}x_p = \sum_{j=1}^p a_{1j}x_j \]
El análisis de componentes principales continúa encontrando una función lineal \(a_2^\prime y\) que no está correlacionada con \(a_1^\prime y\) con máxima varianza y así sucesivamente hasta \(k\) componentes principales.
Ejemplo:
Vamos a realizar un análisis de componentes principales sobre los resultados obtenidos en la competición de heptathlon femenino en los juegos Olímpicos de Seúl (1988).
library(HSAUR2)## Loading required package: tools
data("heptathlon")
head(heptathlon)## hurdles highjump shot run200m longjump javelin
## Joyner-Kersee (USA) 12.69 1.86 15.80 22.56 7.27 45.66
## John (GDR) 12.85 1.80 16.23 23.65 6.71 42.56
## Behmer (GDR) 13.20 1.83 14.20 23.10 6.68 44.54
## Sablovskaite (URS) 13.61 1.80 15.23 23.92 6.25 42.78
## Choubenkova (URS) 13.51 1.74 14.76 23.93 6.32 47.46
## Schulz (GDR) 13.75 1.83 13.50 24.65 6.33 42.82
## run800m score
## Joyner-Kersee (USA) 128.51 7291
## John (GDR) 126.12 6897
## Behmer (GDR) 124.20 6858
## Sablovskaite (URS) 132.24 6540
## Choubenkova (URS) 127.90 6540
## Schulz (GDR) 125.79 6411
Recodificamos las pruebas relativas a 3 carreras hurdles, run200m y run800m, restando al mayor valor en cada carrera, cada uno de los tiempos de las 35 atletas.
heptathlon$hurdles <- max(heptathlon$hurdles) - heptathlon$hurdles
heptathlon$run200m <- max(heptathlon$run200m) - heptathlon$run200m
heptathlon$run800m <- max(heptathlon$run800m) - heptathlon$run800mDiagrama de dispersion
score <- which(colnames(heptathlon) == "score")
plot(heptathlon[,-score])Matriz de correlaciones
round(cor(heptathlon[,-score]),3)## hurdles highjump shot run200m longjump javelin run800m
## hurdles 1.000 0.811 0.651 0.774 0.912 0.008 0.779
## highjump 0.811 1.000 0.441 0.488 0.782 0.002 0.591
## shot 0.651 0.441 1.000 0.683 0.743 0.269 0.420
## run200m 0.774 0.488 0.683 1.000 0.817 0.333 0.617
## longjump 0.912 0.782 0.743 0.817 1.000 0.067 0.700
## javelin 0.008 0.002 0.269 0.333 0.067 1.000 -0.020
## run800m 0.779 0.591 0.420 0.617 0.700 -0.020 1.000
el resultado confirma que la gran mayoría de las correlaciones entre las pruebas son positivas, con una correlación alta entre el salto de longitud (longjump) y los 100m vallas (hurdles). Algunos menos como el salto de altura (highjump) y el lanzamiento de peso (shot) y la jabalina (javelin) que presenta una correlación cercana a cero con el resto de pruebas.
Una posible expliación a este resultado puede deberse a que el entrenamiento para las otras 6 pruebas, no aporta demasiado al de la prueba de jabalina, que es una prueba más técnica.
Se puede observar que hay un valor atípico en casi todas las pruebas que corresponde a una atleta Launa (PNG) de Papua Nueva Guinea- Eliminaremos esta observación para ver si la matriz de correlaciones es sensiblemente distinta:
heptathlon <- heptathlon[-which(rownames(heptathlon)=="Launa (PNG)"),]
plot(heptathlon[,-score])Eliminando a la atleta de Papua Nueva Guinea, las correlaciones cambian sustancialmente y en el diagrama de dispersión de la matriz, no se observan valores extremos.
round(cor(heptathlon[,-score]),3)## hurdles highjump shot run200m longjump javelin run800m
## hurdles 1.000 0.582 0.767 0.830 0.889 0.332 0.559
## highjump 0.582 1.000 0.465 0.391 0.663 0.348 0.152
## shot 0.767 0.465 1.000 0.669 0.784 0.343 0.408
## run200m 0.830 0.391 0.669 1.000 0.811 0.471 0.573
## longjump 0.889 0.663 0.784 0.811 1.000 0.287 0.523
## javelin 0.332 0.348 0.343 0.471 0.287 1.000 0.256
## run800m 0.559 0.152 0.408 0.573 0.523 0.256 1.000
Para realizar el PCA, partiremos de la matriz de correlaciones, ya que las 7 pruebas están medidas en diferentes escalas (metros, segundos). A este procedimiento se le denomina PCA normalizado (scale=TRUE en la función prcomp).
?prcomp
heptathlon_pca <- prcomp(heptathlon[,-score],scale=TRUE)
head(heptathlon_pca)## $sdev
## [1] 2.0793370 0.9481532 0.9109016 0.6831967 0.5461888 0.3374549 0.2620420
##
## $rotation
## PC1 PC2 PC3 PC4 PC5
## hurdles -0.4503876 0.05772161 -0.1739345 0.04840598 -0.19889364
## highjump -0.3145115 -0.65133162 -0.2088272 -0.55694554 0.07076358
## shot -0.4024884 -0.02202088 -0.1534709 0.54826705 0.67166466
## run200m -0.4270860 0.18502783 0.1301287 0.23095946 -0.61781764
## longjump -0.4509639 -0.02492486 -0.2697589 -0.01468275 -0.12151793
## javelin -0.2423079 -0.32572229 0.8806995 0.06024757 0.07874396
## run800m -0.3029068 0.65650503 0.1930020 -0.57418128 0.31880178
## PC6 PC7
## hurdles 0.84665086 -0.06961672
## highjump -0.09007544 0.33155910
## shot -0.09886359 0.22904298
## run200m -0.33279359 0.46971934
## longjump -0.38294411 -0.74940781
## javelin 0.07193437 -0.21108138
## run800m -0.05217664 0.07718616
##
## $center
## hurdles highjump shot run200m longjump javelin run800m
## 2.687500 1.793750 13.173333 2.023750 6.205417 41.278333 28.516667
##
## $scale
## hurdles highjump shot run200m longjump javelin
## 0.51456398 0.05232112 1.49714995 0.93676972 0.40165938 3.46870690
## run800m
## 6.14724800
##
## $x
## PC1 PC2 PC3 PC4
## Joyner-Kersee (USA) -4.757530189 -0.13986143 -0.006040526 0.293416339
## John (GDR) -3.147943402 0.94859029 -0.243919842 0.549171385
## Behmer (GDR) -2.926184760 0.69534239 0.622293440 -0.554744912
## Sablovskaite (URS) -1.288135516 0.17900713 0.250632380 0.637174187
## Choubenkova (URS) -1.503450994 0.96177329 1.780588549 0.784035325
## Schulz (GDR) -0.958467101 0.35121643 0.413086366 -1.113546938
## Fleming (AUS) -0.953445060 0.49982537 -0.265135015 -0.140202490
## Greiner (USA) -0.633239267 0.37592917 -1.140338594 0.142558348
## Lajbnerova (CZE) -0.381571974 -0.71213459 -0.068395353 0.087212735
## Bouraga (URS) -0.522322004 0.77688861 -0.481071429 0.283745698
## Wijnsma (HOL) -0.217701500 -0.23369645 -1.154221444 -1.260128609
## Dimitrova (BUL) -1.075984276 0.51552998 -0.312458252 -0.127032432
## Scheider (SWI) 0.003014986 -1.44688825 1.582739069 -1.254415325
## Braun (FRG) 0.109183759 -1.63595645 0.469577294 0.362580442
## Ruotsalainen (FIN) 0.208868056 -0.68866173 1.152140223 -0.112914470
## Yuping (CHN) 0.232507119 -1.95999641 -1.541230813 0.598325122
## Hagger (GB) 0.659520046 -0.08775813 -1.796509771 -0.182375000
## Brown (USA) 0.756854602 -2.04292201 0.451506018 0.476926314
## Mulliner (GB) 1.880932819 0.91530324 -0.359311801 0.799619094
## Hautenauve (BEL) 1.828170404 0.72629699 -1.048640439 -0.711793161
## Kytola (FIN) 2.118203163 0.39921397 0.190158154 -0.788445056
## Geremias (BRA) 2.770706272 0.03463584 0.170274969 1.385562494
## Hui-Ing (TAI) 3.901166920 1.20175472 0.943677497 -0.002429122
## Jeong-Mi (KOR) 3.896847898 0.36656804 0.390599321 -0.152299968
## PC5 PC6 PC7
## Joyner-Kersee (USA) -0.36183307 -0.27050283 -0.47587527
## John (GDR) 0.75364464 0.37770017 -0.05172711
## Behmer (GDR) -0.19035037 -0.25780287 0.11054960
## Sablovskaite (URS) 0.60362153 -0.21575716 0.53075152
## Choubenkova (URS) 0.58969949 0.08014332 -0.30081842
## Schulz (GDR) 0.71483887 -0.25436956 0.03838796
## Fleming (AUS) -0.86581530 0.03691813 0.23005943
## Greiner (USA) 0.20807431 -0.14236240 -0.06374657
## Lajbnerova (CZE) 0.67727618 0.25014881 0.35555639
## Bouraga (URS) -1.18784299 0.39881271 0.19712215
## Wijnsma (HOL) 0.37497195 -0.20267731 0.17459647
## Dimitrova (BUL) -0.91992929 0.26727067 0.21111846
## Scheider (SWI) -0.20526249 0.17597425 -0.03915701
## Braun (FRG) -0.14712208 0.26134199 -0.01334416
## Ruotsalainen (FIN) -0.31539746 0.18351622 -0.14127555
## Yuping (CHN) 0.17451428 -0.50175724 0.04999374
## Hagger (GB) -0.05104049 0.55058471 -0.46388534
## Brown (USA) -0.38154294 -0.26606429 -0.11099445
## Mulliner (GB) -0.06942955 -0.73259727 -0.31281502
## Hautenauve (BEL) 0.14092347 0.06933542 -0.07548638
## Kytola (FIN) 0.41815113 -0.03363651 0.12143219
## Geremias (BRA) 0.28541366 0.38083979 0.34574480
## Hui-Ing (TAI) -0.67080776 -0.52756760 0.09436975
## Jeong-Mi (KOR) 0.42524426 0.37250885 -0.41055719
a1 <- heptathlon_pca$rotation[,1]Podemos resumir el análisis con summary
summary(heptathlon_pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.0793 0.9482 0.9109 0.68320 0.54619 0.33745
## Proportion of Variance 0.6177 0.1284 0.1185 0.06668 0.04262 0.01627
## Cumulative Proportion 0.6177 0.7461 0.8646 0.93131 0.97392 0.99019
## PC7
## Standard deviation 0.26204
## Proportion of Variance 0.00981
## Cumulative Proportion 1.00000
Los primeros 2 componentes explican un 75% de la varianza explicada y las 3 primeras un 86%.
head(heptathlon_pca$x)## PC1 PC2 PC3 PC4
## Joyner-Kersee (USA) -4.7575302 -0.1398614 -0.006040526 0.2934163
## John (GDR) -3.1479434 0.9485903 -0.243919842 0.5491714
## Behmer (GDR) -2.9261848 0.6953424 0.622293440 -0.5547449
## Sablovskaite (URS) -1.2881355 0.1790071 0.250632380 0.6371742
## Choubenkova (URS) -1.5034510 0.9617733 1.780588549 0.7840353
## Schulz (GDR) -0.9584671 0.3512164 0.413086366 -1.1135469
## PC5 PC6 PC7
## Joyner-Kersee (USA) -0.3618331 -0.27050283 -0.47587527
## John (GDR) 0.7536446 0.37770017 -0.05172711
## Behmer (GDR) -0.1903504 -0.25780287 0.11054960
## Sablovskaite (URS) 0.6036215 -0.21575716 0.53075152
## Choubenkova (URS) 0.5896995 0.08014332 -0.30081842
## Schulz (GDR) 0.7148389 -0.25436956 0.03838796
Graficamente, se observa que la primera componente principal es la más dominante.
plot(heptathlon_pca)El biplot es una representación gráfica de datos multivariantes. De la misma manera que un diagrama de dispersión muestra la distribución conjunta de dos variables, un biplot representa tres o más variables.
Si ordenamos de mayor a menor la variable score tenemos a las tres medallas de oro, plata y bronce.
head(heptathlon[order(heptathlon$score,decreasing = TRUE),],3)## hurdles highjump shot run200m longjump javelin
## Joyner-Kersee (USA) 3.73 1.86 15.80 4.05 7.27 45.66
## John (GDR) 3.57 1.80 16.23 2.96 6.71 42.56
## Behmer (GDR) 3.22 1.83 14.20 3.51 6.68 44.54
## run800m score
## Joyner-Kersee (USA) 34.92 7291
## John (GDR) 37.31 6897
## Behmer (GDR) 39.23 6858
El biplot, nos muestra a las atletas proyectadas sobres sus 2 primeras componentes principales, pero además las flechas nos dan información de las varianzas y covarianzas de las variables (direcciones de máxima variabilidad).
biplot(heptathlon_pca) Por ejemplo, la ganadora de la prueba
Joyner-Kersee (USA) acumula mayores puntuaciones en las pruebas longjump, hurdles y run200m.
Podemos analizar la correlación entre la variable score y la PC1. Lo cual indica que la correlación es muy negativa y muy fuerte con el score.
cor(heptathlon$score, heptathlon_pca$x[,1])## [1] -0.9931168
plot(heptathlon$score, heptathlon_pca$x[,1])Los objetos de la funcion prcomp
class(heptathlon_pca)## [1] "prcomp"
permiten utilizar la funcion predict. Por ejemplo:
new.athlete <- as.data.frame(t(as.vector(c(3.5,2,13,5,7,41,33))))
colnames(new.athlete) <- c("hurdles","highjump",
"shot","run200m","longjump",
"javelin","run800m")
rownames(new.athlete) <- "Mrs XYZ (XXX)"
pp<-predict(heptathlon_pca,newdata = new.athlete)
pp## PC1 PC2 PC3 PC4 PC5 PC6
## Mrs XYZ (XXX) -4.354879 -1.430366 -1.130194 -1.901377 -2.089963 -0.8654821
## PC7
## Mrs XYZ (XXX) 1.253643
biplot(heptathlon_pca)
text(pp[,1],pp[,2],"Mrs XYZ (XXX)",col="blue",cex=.65)Volvamos al ejemplo housing,data. Vamos a tratar de predecir el valor medio de las casas ocupadas por sus dueños en miles de dólares (MEDV) usando las tres primeras PCs de un PCA.
houses <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data",
header = F, na.string = "?")
colnames(houses) <- c("CRIM", "ZN", "INDUS","CHAS",
"NOX","RM","AGE","DIS","RAD",
"TAX","PTRATIO","B","LSTAT","MEDV")
# Perform PCA
pcaHouses <- prcomp(houses[,-14],scale=TRUE)
summary(pcaHouses)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.4752 1.1972 1.11473 0.92605 0.91368 0.81081
## Proportion of Variance 0.4713 0.1103 0.09559 0.06597 0.06422 0.05057
## Cumulative Proportion 0.4713 0.5816 0.67713 0.74310 0.80732 0.85789
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.73168 0.62936 0.5263 0.46930 0.43129 0.41146
## Proportion of Variance 0.04118 0.03047 0.0213 0.01694 0.01431 0.01302
## Cumulative Proportion 0.89907 0.92954 0.9508 0.96778 0.98209 0.99511
## PC13
## Standard deviation 0.25201
## Proportion of Variance 0.00489
## Cumulative Proportion 1.00000
scoresHouses <- pcaHouses$x
# Fit lm using the first 3 PCs
modHouses <- lm(houses$MEDV ~ scoresHouses[,1:3])
summary(modHouses)##
## Call:
## lm(formula = houses$MEDV ~ scoresHouses[, 1:3])
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.774 -3.082 -0.651 2.091 35.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.5328 0.2474 91.06 <2e-16 ***
## scoresHouses[, 1:3]PC1 -2.2730 0.1001 -22.72 <2e-16 ***
## scoresHouses[, 1:3]PC2 2.1949 0.2069 10.61 <2e-16 ***
## scoresHouses[, 1:3]PC3 3.5010 0.2222 15.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.566 on 502 degrees of freedom
## Multiple R-squared: 0.6359, Adjusted R-squared: 0.6337
## F-statistic: 292.3 on 3 and 502 DF, p-value: < 2.2e-16
Comparamos con el modelo lineal con todas las variables
# Fit lm using all 14 vars
modHousesFull <- lm(MEDV ~ ., data = houses)
summary(modHousesFull) # R2 = 0.741##
## Call:
## lm(formula = MEDV ~ ., data = houses)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.595 -2.730 -0.518 1.777 26.199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.646e+01 5.103e+00 7.144 3.28e-12 ***
## CRIM -1.080e-01 3.286e-02 -3.287 0.001087 **
## ZN 4.642e-02 1.373e-02 3.382 0.000778 ***
## INDUS 2.056e-02 6.150e-02 0.334 0.738288
## CHAS 2.687e+00 8.616e-01 3.118 0.001925 **
## NOX -1.777e+01 3.820e+00 -4.651 4.25e-06 ***
## RM 3.810e+00 4.179e-01 9.116 < 2e-16 ***
## AGE 6.922e-04 1.321e-02 0.052 0.958229
## DIS -1.476e+00 1.995e-01 -7.398 6.01e-13 ***
## RAD 3.060e-01 6.635e-02 4.613 5.07e-06 ***
## TAX -1.233e-02 3.760e-03 -3.280 0.001112 **
## PTRATIO -9.527e-01 1.308e-01 -7.283 1.31e-12 ***
## B 9.312e-03 2.686e-03 3.467 0.000573 ***
## LSTAT -5.248e-01 5.072e-02 -10.347 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
## F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16
# Compare obs. vs. pred. plots
par(mfrow = c(1,2))
plot(houses$MEDV, predict(modHouses),
xlab = "Observed MEDV", ylab = "Predicted MEDV", main = "PCR", abline(a = 0, b = 1, col = "red"))
plot(houses$MEDV, predict(modHousesFull),
xlab = "Observed MEDV", ylab = "Predicted MEDV", main = "Full model", abline(a = 0, b = 1, col = "red"))Ejemplo
data(USArrests)
head(USArrests)## Murder Assault UrbanPop Rape
## Alabama 13.2 236 58 21.2
## Alaska 10.0 263 48 44.5
## Arizona 8.1 294 80 31.0
## Arkansas 8.8 190 50 19.5
## California 9.0 276 91 40.6
## Colorado 7.9 204 78 38.7
# PCA with function prcomp
pca1 = prcomp(USArrests, scale. = TRUE)
# sqrt of eigenvalues
pca1$sdev## [1] 1.5748783 0.9948694 0.5971291 0.4164494
# loadings
head(pca1$rotation)## PC1 PC2 PC3 PC4
## Murder -0.5358995 0.4181809 -0.3412327 0.64922780
## Assault -0.5831836 0.1879856 -0.2681484 -0.74340748
## UrbanPop -0.2781909 -0.8728062 -0.3780158 0.13387773
## Rape -0.5434321 -0.1673186 0.8177779 0.08902432
# PCs (aka scores)
head(pca1$x)## PC1 PC2 PC3 PC4
## Alabama -0.9756604 1.1220012 -0.43980366 0.154696581
## Alaska -1.9305379 1.0624269 2.01950027 -0.434175454
## Arizona -1.7454429 -0.7384595 0.05423025 -0.826264240
## Arkansas 0.1399989 1.1085423 0.11342217 -0.180973554
## California -2.4986128 -1.5274267 0.59254100 -0.338559240
## Colorado -1.4993407 -0.9776297 1.08400162 0.001450164
plot(pca1)biplot(pca1)Podemos cambiar el eje
pca.out <- pca1
pca.out$rotation <- -pca.out$rotation
pca.out$x <- -pca.out$x
biplot(pca.out,scale=0, cex=.7)Interpretación del biplot: Para cada una de los 50 estados en los Estados Unidos, el conjunto de datos contiene el número de arrestos por cada 100.000 residentes por cada tres crímenes: Asalto, Asesinato y Violación. También Urbanpop representa el porcentaje de la población en cada estado que vive en áreas urbanas. La gráfica muestra las dos primeras puntuaciones de los componentes principales y el vector loadings en gráfico simple.
pca.out$rotation[,1:2]## PC1 PC2
## Murder 0.5358995 -0.4181809
## Assault 0.5831836 -0.1879856
## UrbanPop 0.2781909 0.8728062
## Rape 0.5434321 0.1673186
El PC1, muestra como el peso de las variables Murder, Assault y Rape es similar, y más pequeño para UrbanPop. Por lo tanto, este componente corresponde aproximadamente a una medida de las tasas globales de delitos graves.
El PC2, tiene un peso mucho mayor en Urbanpop, y mucho menor en el resto. Por lo tanto, este componente corresponde aproximadamente al nivel de urbanización del estado.
En general, vemos que las variables relacionadas con el delito están situadas cerca unas de otras, y que la variable urbanpop está lejos de las otras tres. Esto indica que las variables relacionadas con el delito están correlacionadas entre sí. Los estados con altas tasas de homicidio tienden a tener altas tasas de agresión y violación. La variable Urbanpop está menos correlacionada con las otras tres. Otro paquete de R es FactorMineR
# PCA with function PCA
library(FactoMineR)## Warning: package 'FactoMineR' was built under R version 3.3.2
# apply PCA
pca3 = PCA(USArrests, graph = FALSE)
# matrix with eigenvalues
pca3$eig## eigenvalue percentage of variance cumulative percentage of variance
## comp 1 2.4802416 62.006039 62.00604
## comp 2 0.9897652 24.744129 86.75017
## comp 3 0.3565632 8.914080 95.66425
## comp 4 0.1734301 4.335752 100.00000
# correlations between variables and PCs
pca3$var$coord## Dim.1 Dim.2 Dim.3 Dim.4
## Murder 0.8439764 -0.4160354 0.2037600 0.27037052
## Assault 0.9184432 -0.1870211 0.1601192 -0.30959159
## UrbanPop 0.4381168 0.8683282 0.2257242 0.05575330
## Rape 0.8558394 0.1664602 -0.4883190 0.03707412
plot(pca3)Ejemplo
data("USairpollution")
head(USairpollution)## SO2 temp manu popul wind precip predays
## Albany 46 47.6 44 116 8.8 33.36 135
## Albuquerque 11 56.8 46 244 8.9 7.77 58
## Atlanta 24 61.5 368 497 9.1 48.34 115
## Baltimore 47 55.0 625 905 9.6 41.31 111
## Buffalo 11 47.1 391 463 12.4 36.11 166
## Charleston 31 55.2 35 71 6.5 40.75 148
library(TeachingDemos)
faces(USairpollution)faces2(USairpollution)panel.hist <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col="blue", ...)
}
pairs(USairpollution,diag.panel=panel.hist)# Calculo la matriz de correlaciones
cor(USairpollution[,-1]) ## temp manu popul wind precip
## temp 1.00000000 -0.19004216 -0.06267813 -0.34973963 0.38625342
## manu -0.19004216 1.00000000 0.95526935 0.23794683 -0.03241688
## popul -0.06267813 0.95526935 1.00000000 0.21264375 -0.02611873
## wind -0.34973963 0.23794683 0.21264375 1.00000000 -0.01299438
## precip 0.38625342 -0.03241688 -0.02611873 -0.01299438 1.00000000
## predays -0.43024212 0.13182930 0.04208319 0.16410559 0.49609671
## predays
## temp -0.43024212
## manu 0.13182930
## popul 0.04208319
## wind 0.16410559
## precip 0.49609671
## predays 1.00000000
Otra funcion para calcular PCA en princomp
# Calculo los componentes principales basados en la matriz de correlaciones
USair.pc<-princomp(USairpollution[,-1],cor=TRUE)
summary(USair.pc,loadings=TRUE)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 1.4819456 1.2247218 1.1809526 0.8719099 0.33848287
## Proportion of Variance 0.3660271 0.2499906 0.2324415 0.1267045 0.01909511
## Cumulative Proportion 0.3660271 0.6160177 0.8484592 0.9751637 0.99425879
## Comp.6
## Standard deviation 0.185599752
## Proportion of Variance 0.005741211
## Cumulative Proportion 1.000000000
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## temp 0.330 -0.128 0.672 -0.306 0.558 -0.136
## manu -0.612 -0.168 0.273 0.137 -0.102 -0.703
## popul -0.578 -0.222 0.350 0.695
## wind -0.354 0.131 -0.297 -0.869 0.113
## precip 0.623 0.505 -0.171 -0.568
## predays -0.238 0.708 0.311 0.580
#scores
head(USair.pc$scores[,1:3])## Comp.1 Comp.2 Comp.3
## Albany 0.5389456 0.79206938 -1.3626052
## Albuquerque 1.4170931 -2.86589024 -1.2754003
## Atlanta 0.5989948 0.58723563 0.9954113
## Baltimore -0.5093800 0.02875397 0.3635995
## Buffalo -1.3907089 1.88030104 -1.7761919
## Charleston 1.4297533 1.21058365 0.0794435
# PC1 vs PC2
par(pty="s")
plot(USair.pc$scores[,1],USair.pc$scores[,2],
ylim=range(USair.pc$scores[,1]),
xlab="PC1",ylab="PC2",type="n",lwd=2)
text(USair.pc$scores[,1],USair.pc$scores[,2],
labels=abbreviate(row.names(USairpollution)),cex=0.7,lwd=2)# PC1 vs PC3
par(pty="s")
plot(USair.pc$scores[,1],USair.pc$scores[,3],
ylim=range(USair.pc$scores[,1]),
xlab="PC1",ylab="PC3",type="n",lwd=2)
text(USair.pc$scores[,1],USair.pc$scores[,3],
labels=abbreviate(row.names(USairpollution)),cex=0.7,lwd=2)# PC2 vs PC3
par(pty="s")
plot(USair.pc$scores[,2],USair.pc$scores[,3],
ylim=range(USair.pc$scores[,2]),
xlab="PC1",ylab="PC3",type="n",lwd=2)
text(USair.pc$scores[,2],USair.pc$scores[,3],
labels=abbreviate(row.names(USairpollution)),cex=0.7,lwd=2)El análisis de cluster es una técnica cuya idea básica es agrupar un conjunto de observaciones en un número dado de clusters o grupos. Este agrupamiento se basa en la idea de distancia o similitud entre las observaciones. Por tanto, la obtención de dichos clusters depende del criterio o distancia considerados.
Ejemplos de distancias:
\[ d_{I_1,I_2} = (x_i - x_j)^{\prime} W^{-1} (x_i - x_j) \] donde \(W\) es la matriz de covarianzas entre las variables. De este modo, las variables se ponderan según el grado de relación que exista entre ellas, es decir, si están más o menos correlacionadas. Si la correlación es nula y las variables están estandarizadas, se obtiene la distancia euclídea.
En la práctica, no se pueden examinar todas las posibilidades de agrupar los elementos, incluso con los ordenadores más rápidos. Una solución se encuentra en los llamados métodos jerárquicos.
Métodos aglomerativos: se comienza con los objetos o individuos de modo individual; de este modo, se tienen tantos clusters iniciales como objetos. Luego se van agrupando de modo que los primeros en hacerlo son los más similares y al final, todos los subgrupos se unen en un único cluster.
Métodos divisivos o de particionado: al contrario. Se parte de un grupo único con todas las observaciones y se van dividiendo según lo lejanos que estén.
En cualquier caso, de ambos métodos se deriva un dendograma, que es un gráfico que ilustra cómo se van haciendo las subdivisiones o los agrupamientos, etapa a etapa.
Consideramos aquí los métodos aglomerativos con diferentes métodos de unión (linkage methods). Los más importantes son:
Mínima distancia o vecino más próximo.
Máxima distancia o vecino más lejano.
Distancia media (average distance).
Se puede observar que, de este modo, se define una posible distancia entre dos clusters: la correspondiente a la pareja de elementos más cercana, la más lejana o la media de todas las posibles parejas de elementos de ambos clusters.
Se usan para agrupar objetos, pero no variables, en un conjunto de \(k\) clusters ya predeterminado. No se tiene que especificar una matriz de distancias, lo cual entre otras razones permite trabajar con un número de datos mayor que en el caso de los métodos jerárquicos.
Se parte de un conjunto inicial de clusters elegidos al azar, que son los representantes de todos ellos; luego se van cambiando de modo iterativo. Se usa habitualmente el método de las \(k\)-medias.
Método de las k-medias permite asignar a cada observación el cluster que se encuentra más próximo en términos del centroide (media). En general, la distancia empleada es la euclídea.
Pasos:
Se toman al azar \(k\) clusters iniciales.
Para el conjunto de observaciones, se vuelve a calcular las distancias a los centroides de los clusters y se reasignan a los que estén más próximos. Se vuelven a recalcular los centroides de los \(k\) clusters después de las reasignaciones de los elementos.
Se repiten los dos pasos anteriores hasta que no se produzca ninguna reasignación, es decir, hasta que los elementos se estabilicen en algún grupo. Usualmente, se especifican \(k\) centroides iniciales y se procede al paso (2) y, en la práctica, se observan la mayor parte de reasignaciones en las primeras iteraciones.
### Ejemplos:
En este conjunto de datos se observa la composición de diferentes vinos.
# install.packages('rattle')
data(wine, package='rattle')
head(wine)## Type Alcohol Malic Ash Alcalinity Magnesium Phenols Flavanoids
## 1 1 14.23 1.71 2.43 15.6 127 2.80 3.06
## 2 1 13.20 1.78 2.14 11.2 100 2.65 2.76
## 3 1 13.16 2.36 2.67 18.6 101 2.80 3.24
## 4 1 14.37 1.95 2.50 16.8 113 3.85 3.49
## 5 1 13.24 2.59 2.87 21.0 118 2.80 2.69
## 6 1 14.20 1.76 2.45 15.2 112 3.27 3.39
## Nonflavanoids Proanthocyanins Color Hue Dilution Proline
## 1 0.28 2.29 5.64 1.04 3.92 1065
## 2 0.26 1.28 4.38 1.05 3.40 1050
## 3 0.30 2.81 5.68 1.03 3.17 1185
## 4 0.24 2.18 7.80 0.86 3.45 1480
## 5 0.39 1.82 4.32 1.04 2.93 735
## 6 0.34 1.97 6.75 1.05 2.85 1450
Siempre es recomendable estandarizar las variables
wine.stand <- scale(wine[-1]) # To standarize the variables
# K-Means
k.means.fit <- kmeans(wine.stand, 3) # k = 3attributes(k.means.fit)## $names
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
##
## $class
## [1] "kmeans"
## Alcohol Malic Ash Alcalinity Magnesium Phenols
## 1 0.1644436 0.8690954 0.1863726 0.5228924 -0.07526047 -0.97657548
## 2 -0.9234669 -0.3929331 -0.4931257 0.1701220 -0.49032869 -0.07576891
## 3 0.8328826 -0.3029551 0.3636801 -0.6084749 0.57596208 0.88274724
## Flavanoids Nonflavanoids Proanthocyanins Color Hue
## 1 -1.21182921 0.72402116 -0.77751312 0.9388902 -1.1615122
## 2 0.02075402 -0.03343924 0.05810161 -0.8993770 0.4605046
## 3 0.97506900 -0.56050853 0.57865427 0.1705823 0.4726504
## Dilution Proline
## 1 -1.2887761 -0.4059428
## 2 0.2700025 -0.7517257
## 3 0.7770551 1.1220202
## [1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [36] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 1 2 2 2 2 2 2 2 2
## [71] 2 2 2 3 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2
## [106] 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 3 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1
## [141] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [176] 1 1 1
## [1] 51 65 62
Una pregunta fundamental es cómo determinar el valor del parámetro \(k\). Si consideramos el porcentaje de varianza explicado como una función del número de grupos:
Uno debe elegir un número de grupos de modo que la adición de otro grupo no da mucho mejor modelado de los datos. Más precisamente, si se traza el porcentaje de varianza explicado por los conglomerados en función del número de racimos, los primeros grupos agregarán mucha información (explican mucha varianza), pero en algún momento la ganancia marginal disminuirá, dando un ángulo en la gráfico. El número de racimos se elige en este punto, de ahí el “criterio del codo”.
wssplot <- function(data, nc=15, seed=1234){
wss <- (nrow(data)-1)*sum(apply(data,2,var))
for (i in 2:nc){
set.seed(seed)
wss[i] <- sum(kmeans(data, centers=i)$withinss)}
plot(1:nc, wss, type="b", xlab="Number of Clusters",
ylab="Within groups sum of squares")}
wssplot(wine.stand, nc=6) La library(cluster) permite (con el paso previo de un PCA) representar los clusters en 2 dimensiones
library(cluster)
clusplot(wine.stand, k.means.fit$cluster,
main='2D representation of the Cluster solution',
color=TRUE, shade=TRUE,
labels=2, lines=0)Sabiendo que hay tres tipos de vinos wine$Type, podemos calcular una matriz de confusión
table(wine$Type)##
## 1 2 3
## 59 71 48
table(wine[,1],k.means.fit$cluster)##
## 1 2 3
## 1 0 0 59
## 2 3 65 3
## 3 48 0 0
Cluster jerarquico
Los métodos jerárquicos utilizan una matriz de distancia como una input para el algoritmo de agrupación. La elección de una métrica apropiada influirá en la forma de los clusters, ya que algunos elementos pueden estar cerca uno del otro según una distancia y más lejos de acuerdo con otro.
d <- dist(wine.stand, method = "euclidean") # Euclidean distance matrix.Ward’s minimum variance criterion minimizes the total within-cluster variance
H.fit <- hclust(d, method="ward")## The "ward" method has been renamed to "ward.D"; note new "ward.D2"
class(H.fit)## [1] "hclust"
La opcion plot de un objeto hclust nos devuelve el dendograma:
plot(H.fit) # display dendogram
groups <- cutree(H.fit, k=3) # cut tree into 5 clusters
# draw dendogram with red borders around the 5 clusters
rect.hclust(H.fit, k=3, border="red") table(wine[,1],groups)## groups
## 1 2 3
## 1 58 1 0
## 2 7 58 6
## 3 0 0 48
R tiene amplias facilidades para analizar datos de series temporales. Esta sección describe la creación de una serie de tiempo, descomposición estacional, predicción con el paquete forecast.
# save a numeric vector containing 84 monthly observations
# from Jan 2009 to Dec 2015 as a time series object
set.seed(1234)
myvector <- rnorm(84)
myts <- ts(myvector, start=c(2009, 1), end=c(2015, 12), frequency=12)
class(myts)## [1] "ts"
# subset the time series (June 2014 to December 2014)
myts2 <- window(myts, start=c(2014, 6), end=c(2014, 12))
# plot series
plot(myts)Ejemplo con datos reales:
data(AirPassengers)
class(AirPassengers)## [1] "ts"
?AirPassengers
start(AirPassengers)## [1] 1949 1
end(AirPassengers)## [1] 1960 12
frequency(AirPassengers)## [1] 12
cycle(AirPassengers)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 1 2 3 4 5 6 7 8 9 10 11 12
## 1950 1 2 3 4 5 6 7 8 9 10 11 12
## 1951 1 2 3 4 5 6 7 8 9 10 11 12
## 1952 1 2 3 4 5 6 7 8 9 10 11 12
## 1953 1 2 3 4 5 6 7 8 9 10 11 12
## 1954 1 2 3 4 5 6 7 8 9 10 11 12
## 1955 1 2 3 4 5 6 7 8 9 10 11 12
## 1956 1 2 3 4 5 6 7 8 9 10 11 12
## 1957 1 2 3 4 5 6 7 8 9 10 11 12
## 1958 1 2 3 4 5 6 7 8 9 10 11 12
## 1959 1 2 3 4 5 6 7 8 9 10 11 12
## 1960 1 2 3 4 5 6 7 8 9 10 11 12
plot(AirPassengers)¿Qué patrones ves?
La tendencia interanual muestra claramente que el número de pasajeros ha ido en aumento.
# plot of aggregated by months
plot(aggregate(AirPassengers,FUN=mean)) La varianza y el valor medio en julio y agosto es mucho mayor que el resto de los meses.
# boxplots accross months will show seasonal effects
boxplot(AirPassengers~cycle(AirPassengers)) También observamos un fuerte efecto estacional con un ciclo de 12 meses o menos.
La función decompose() descompone una serie de tiempo en la estacionalidad, la tendencia y los componentes irregulares utilizando promedios móviles. El componente estacional puede ser aditivo: \[
Y_t = T_t + S_t + e_t
\]
forecastlibrary(forecast)
tsData <- EuStockMarkets[, 1] # ts data
decomposedRes <- decompose(tsData,
type="mult")
# use type = "additive" for additive components
plot (decomposedRes) # see plot belowstlRes <- stl(tsData, s.window = "periodic")library(forecast)
ts.stl <- stl(AirPassengers,"periodic") # decompose the TS
ts.sa <- seasadj(ts.stl) # de-seasonalize
plot(AirPassengers, type="l") # original seriesplot(ts.sa, type="l") # seasonal adjustedseasonplot(ts.sa, 12, col=rainbow(12), year.labels=TRUE,
main="Seasonal plot: Airpassengers") # seasonal frequency set as 12 for monthly data.library(forecast)
plot(forecast(ts.stl,h=10))